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Abstract 


A transfer  matrix  for  CO 


laser  beams  with  an  assvimed 


gaussian  intensity  distribution  is  developed  that  Includes 

parameters  for  absorption,  turbulence,  and  thermal  bloomlnB: 

2 < 

Calculated  parameters  are  an  effective  beam  radius  { 1/e 
point)  and  an  on-axls  intensity.  For  moderate  power  levels 


results  are  consistent  with  the  computer  code  COMBO.  The 


blooming  model  predicts  results  worse  than  those  predicted 


by  COMBO  for  high  power  levels 


A MATRIX  APPROACH  TO  A PROPAGATION  CODE 


I.  Introduction 

The  primary  concern  In  nearly  all  high  energy  laser 
applications  Is  the  fraction  of  transmitted  power  that  reaches 
the  receiver.  Calculating  this  fraction  for  atmospheric 
transmission  Is  difficult  due  to  the  Inhomogenlty  of  the  at- 
mosphere as  a medium,  and  to  the  modifications  caused  by  the 
passage  of  the  laser  beam  Itself,  In  addition,  various 
transmitter/receiver  scenarios,  such  as  ground  to  air,  air 
to  air,  stationary  to  moving,  etc.,  will  affect  the  received 
power. 

Thus  there  exists  a need  for  a method  of  easily  deter- 
mining how  efficiently  a given  laser  beam  will  pass  through 
the  atmosphere. 

Background 

There  are  In  existence  several  computer  codes  that  will 
calculate  how  well  a laser  beam  will  propagate  through  the 
atmosphere.  However,  nearly  all  these  codes  are  long  and 
CTmbersome,  requiring  large  amounts  of  computer  storage  and 
computer  time,  or  are  based  on  data  tables  derived  from 
experimental  observations,  again  requiring  large  storage 
capability. 

Both  code  types  require  a computer  to  operate,  are  not 
readily  available,  and  are  not  easily  portable.  Also, 


i 


fe ' 

i, 
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tabular  codes,  while  perhaps  useful,  lack  an  apparent  basis 
In  physics. 

Problem 

The  purpose  of  this  research  is  to  develop  a fast,  easy 
to  apply  method  of  calculating  laser  transmission  that  does 
not  necessarily  require  a computer.  This  method,  or  code, 

I 

Is  to  have  a basis  In  physics  and  Is  designed  to  determine  | 

the  on-axls  Intensity  and  radius  of  the  laser  bieam  at  the  \ 

receiver  when  given  an  Initial  transmitted  power  and 
aperture  radliis. 

Scope  and  Assumptions 

The  following  code  Is  designed  primarily  for  continuous 
wave  CO2  (10,6  lam)  laser  beams.  Only  continuous  wave  beams 
ai*e  considered  so  that  the  beam  can  be  assximed  to  be  In 
steady  state  conditions.  The  areas  of  the  propagation  pro- 
blem Included  here  are  absorption,  turbulence,  and  thermal 
blooming.  A linear  scattering  term  Is  included  In  the 
absorption  section.  The  thermal  blooming  expression  has  a 
correction  factor  for  transmitter  rotation  (beam  slewing) 
and  wind. 

There  are  three  primary  assumptions  made.  The  laser 
beam  Is  assumed  to  have  a gausslan  Intensity  distribution 
at  all  times.  The  atmosphere  Is  modeled  as  a series  of 
homogeneous  layers,  each  layer  having  a thickness  small 
compared  with  the  total  path  length.  Finally,  changes  in 
beam  direction  are  assumed  to  be  small  so  that  paraxial 
approximations  and  geometrical  optics  are  valid. 
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Approach 

An  approach  based  on  ray  transfer  matrices  Is  used.  In 
geometrical  optics  a light  ray  can  be  traced  through  a sys- 
tem using  a transfer  matrix  derived  from  the  system  charac- 
teristics (Ref  4i171-175).  Used  for  yfears  In  designing 
lenses  and  lens  systems,  matrix  techniques  have  not  been 
applied  previously  to  laser  transmission. 

In  this  paper  a separate  matrix  Is  developed  for  each 
part  of  the  propagation  problem.  The  three  Individual 
matrices  are  then  linearly  combined  Into  one  overall  transfer 
matrix.  This  final  matrix  Is  used  to  transfer  the  laser  beam 
through  the  atmospheric  layers  having  predetermined  parameters 
(absorption  coefficients,  wind  velocities,  altitude,  width). 

There  are  several  advantages  to  this  approach.  Matrices 
are  easy  to  use,  especially  In  repeated  applications  of  the 
same  equations.  If  one  area  of  the  propagation  problem  Is 
of  no  Interest,  the  matrix  for  that  area  may  be  deleted  or 
modified  without  affecting  the  other  matrices.  In  addition, 
beam  Intensity  and  size  can  be  determined  at  any  point  along 
the  propagation  path. 

To  determine  If  calculations  using  the  matrix  code  are 
reasonable,  matrix  results  are  compared  with  results  from 
the  computer  code  COMBO  (Ref  6).  COMBO  appears  to  be  an 
aggregate  code  based  on  derived  equations  and  also  tables 
and  curves  formed  from  experimental  data  and  observations. 


since  0 does  not  change  for  this  particular  system.  If  a 
matrix  format  Is  used,  the  equations  can  be  written  as 


I 
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(4) 

or,  In  a similar  way,  equations  can  be  derived  for  any  sys- 
tem or  medltun  and  represented  In  the  form 


""  “ 

- - 

^2 

C D 

0 

E F 

which  Is  the  standard  ray  transfer  matrix  for  geometrical 
optics.  In  this  manner  the  off-axis  position  (radius)  and 
divergence  of  a beam  can  be  determined  for  any  system,  if  6 
Is  small.  Since  lasers  have  beam  divergences  on  the  order 
of  mllllradlans,  this  assumption  Is  valid. 

However,  the  standard  matrix  contains  no  provision  for 
determining  beam  Intensity,  or  equivalently,  beam  power.  To 
represent  power  loss.  It  Is  necessary  to  use  a matrix  of  the 
form 


^2' 

r 1 

ABC 

h 

S 

D E F 

®2 

G H I 

«i 

(6) 


where  the  matrix  elements  A through  I must  now  be  determined. 


Absorption  Matrix 

Linear  absorption  losses  for  CO2  laser  beams  are  primarily 
due  to  aerosol  absorption  and  scattering,  water  absorption, 
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and  CO2  absorption.  For  a given  path  length  between  Initial 
and  final  positions  of  (N  - 1)  and  N,  respectively,  the 
final  beam  power  Is  calculated  from 
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ii 
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i 
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f I 

L ' 

I 
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I 

-I 


(?) 


where  Is  the  absorption  due  to  aerosols,  Is  the  absorp- 
tion due  to  water  vapor,  and  a Is  the  absorption  due  to  at- 

c 

mospherlc  COg.  The  units  of  a are  usually  km“^. 

The  effects  of  aerosol  scattering  and  absorption  are 

strongly  dependent  on  local  atmospheric  conditions.  Location, 

particle  size,  particle  density,  and  other  variables  contribute 

to  aerosol  loss.  See,  for  example,  W.  A.  Proctor's  article, 

"A  Short  Descriptive  Survey  of  Atmospheric  Aerosols"  (Laser 

Digest.  June  I963).  A general  model  Is  beyond  the  scope  of 

this  paper,  but  average  sea  level  coefficients  of  a = O.OI5 

s 

km“^  and  = 0.073  km“^  for  clear  (23  km  visibility)  and 
hazy  (5  km  visibility)  conditions,  respectively,  may  be  used 
(Ref  2 1 1486).  For  high  altitude  paths  or  visibility  greater 
than  25  km,  a « 0 Is  a reasonable  approximation. 

0 

For  later  calculations.  It  Is  convenient  to  define  an 
aerosol  coefficient  b(z)  such  that 


b(z)  s exp(-agZ)  (8) 

where  b(z)  Is  unitless. 

The  water  vapor  coefficient  Is  also  calculated  In 
units  of  km”^.  For  horizontal  beam  paths,  as  a function 
of  altitude  H may  be  calculated  from 
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TABLE  I 


VALUES  OP  A AND  B 


Humidity 

Summer 
(B  = 0.^05) 

(T  t=  30°C) 

Winter 
(B  = 0.575) 

(T  s=  250c) 

A 

A 

{%) 

(mm  Hg) 

(mm  Hg) 

10 

3.28 

2.43 

20 

6.56 

4.85 

30 

9.85 

7.28 

40 

13.13 

9.70 

50 

16.41 

12.12 

60 

19.69 

14.55 

70 

23.78 

16.98 

80 

26.26 

19.40 

90 

29.5^ 

21.83 

a (H)  * 8.3^xlO"^(3.9^AexpC-H(B  + 0.13)]  + A^expC-2BH])  (9) 

where  H Is  In  km  (Ref  5* 1^7^- 1^76) . The  constants  In  equation 
(9)  are  derived  from  experimental  measurements  made  by  the 
authors  of  Ref  5.  The  value  of  A Is  the  water  vapor  pressure 
at  ground  level  In  mm  Hg  and  Is  a function  of  temperature  and 
hvimldlty,  B Is  the  exponential  coefficient  of  the  change  In 
water  vapor  pressure  with  altitude  and  has  the  units  of  km“^ 

B Is  assumed  to  be  a constant  for  summer  and  winter.  Some 
values  for  A and  B are  listed  In  Table  I.  See  Appendix  A 
for  further  explanation. 

For  beam  paths  that  change  altitude,  equation  (9)  Is 
Integrated  over  the  change  In  altitude  to  give  as  a func- 
tion of  final  altitude  and  Initial  altitude  H^.  The 
result  Is 

—4  2 

■■"(expC-2BH^]  - expC-2BH2])  + | 

^|p~^yY3^(expC-Hi(B  + 0.13)]  - expC-HgCB  + 0.13)])  (10) 

where  still  has  the  units  of  km”^  (Ref  5*1476). 

For  atmospheric  CO,  absorption,  the  absorption  coeffl- 

'i 

dents  determined  by  Yin  and  Long  are  used  (Ref  9).  Values  ; 

of  a as  a function  of  altitude  may  be  found  from  the  mraph 
c 

In  Fig.  2.  Appendix  B contains  a list  of  polynomials  used  ] 

to  plot  Fig.  2.  For  beam  paths  that  change  altitude,  the  i 

polynomials  listed  In  Appendix  B are  Integrated  over  the  | 

change  In  altitude  to  give  a as  a function  of  (H-,H, ). 

c cl 

Table  II  lists  some  values  of  a (H-,H,). 

0 2 1 

The  values  of  a (H,,H. ) and  a (H_,H. ) as  previously 
H c 1 C 2 1 
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TABLE  II 


COg  ABSORPTION  COEFFICIENTS 


1 «2 

January 

July 

a 

a 

c 

c 

(km) 

(km"^) 

(km“^) 

0-1 

0.060 

0.079 

0-2 

0.115 

0.148 

0-3 

0.165 

0.208 

0-4 

0.207 

0.261 

0-5 

0.242 

0.307 

0-6 

0.272 

0.347 

0-7 

0.297 

0.381 

0-8 

0.318 

0.409 

0-9 

0.335 

0.433 

0-10 

0.348 

0.451 

0-12 

0.367 

0.477 

0-20 

0.408 

0.520 

(From  Ref  9 i 1553) 
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defined  are  valid  only  for  vertical  paths.  For  slant  paths, 
the  equation  for  the  total  absorption  coefficient  a In  km“^ 


a = csclCa^CHg.H^)  + a^(H2,H^)]  (13 

where  0 Is  the  elevation  angle  with  respect  to  the  horizon 
of  the  laser  beam.  See  Pig.  3. 

Thus  the  final  expression 
for  absorption,  Including 
aerosol  scattering.  Is 


Pjj  = b(z)Pj^_j^exp(-az)  (12) 


where  z la  In  km.  The  Fig.  3.  Beam  Traversing  a 

Slant  Path 

transfer  matrix  for  absorption  only  Is 


Ib(z)exp(-az)  0 


Turbulence 

Atmospheric  turbulence  is  modeled  as  a series  of  homo- 
geneous bubbles,  or  turbules,  of  varying  refractive  Indices. 
As  the  beam  passes  throiU5h  the  turbules.  It  Is  refracted, 
and  the  net  effect  Is  to  Increase  the  beam  radius. 


The  Increase  In  radius  due  to  turbulence  may  be  calculated 


from 


q » 


9 


c/2- 


where  q,  z,  X,  and  p have  the  units  of  m.  The  value  for  p 
Is  calculated  from  the  equation 


P = [57.24  X"^sec0/ (1-^)  dl]"^^^  (I5) 

2 

where  C„  is  the  refractive  Index  structure  constant,  s and  t 
N 

are  the  Initial  and  final  path  points,  and  L = t - s Is  the 

total  length  of  the  propagation  path  In  m (Ref  10i2771). 

2 

The  refractive  Index  structure  constant  C„  Is  a measure  of 

N 

how  the  refractive  Index  varies  due  to  turbulence  as  a 
function  of  altitude.  Several  plots  of  are  shown  In 
Fig.  4. 

2 

Since  Is  not  a function  of  path  length,  equation  (I5) 
can  be  Integrated  to  give 

p = (21.47  X’^LC^  sec0)"^^^  (16) 

2 

An  approximate  expression  for  as  a f\uictlon  of  h Is 

C^(h)  = (17) 

where  h Is  the  mean  altitude  In  m of  the  beam  path  (Ref  61 
20).  This  Is  equivalent  to  letting  h be  represented  by 


h = ^ + ^ (18) 

where  h^  Is  the  altitude  of  the  transmitter  and  Is  set  equal 
to  1 m for  ground  level.  Combining  equations  (14)  through 
(18)  gives  the  expression 
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i 


1 

(19) 


Equation  (19)  Is  Independent  of  the  matrix  parameters,  but  a 
division  by  the  beam  radius  a allows  the  transfer  matrix  for 
Increases  In  the  radius  only  to  be  written  as 


Pxt' 

0 

0 

o" 

• • 

P.T  4 

N 

N-1 

0 

1 + — 
a 

e 

*N-1 

(20)  ! 

0 

0 

0 

j 

N 

Thermal  Bloomlnp; 

Thermal  blooming  Is  a non-llnear  effect  associated  with 
high  energy  laser  beams.  For  CO^  laser  beams  It  IS  mainly 
a result  of  power  absorption  by  atmospheric  CO^  and  water 
vapor.  The  absorbed  power  heats  the  atmosphere  and  causes 
the  air  density  to  decrease  at  the  beam  center.  As  a result 
the  atmosphere  acts  like  a negative  lens,  spreading  the  beam 
and  reducing  on-axls  Intensity. 

The  steady-state  energy  equation  describing  the  power 
absorption  Is 

dCpV  ^ - kV  T = al  (21) 

where  d Is  the  gas  mass  density  In  g/m  , c Is  the  specific 

P 

heat  at  constant  pressuire  In  J/g  °K,  v Is  the  wind  velocity 
In  m/sec  In  the  X direction,  k Is  the  thermal  conductivity 
In  J/m-sec  °K,  T Is  the  temperature  in  °K,  and  I Is  the  Inten- 
slty  In  w/m  (Ref  8i601).  For  no  wind  (v  = 0),  the  develop- 
ment outlined  by  Smith  (Ref  8)  Is  followed.  Letting  I » 


2.00  s. 


L seel?  X 


Lsln0  .1.075 
9 f 


3/5 
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mmmmBm 


P 2 

I^exp(-2r'^/a  ) and  using  the  relation 


d0  l/Snx  /STx 

Fz  “ n'dT^'dr^ 


(22) 


the  angular  divergence  due  to  blooming  Is 


e . e 
r o 


P(an/aT)(l  - expC-2r^/a^1)  (1  - expC-azD 

2Tfknr 


where  0^  is  the  beam  divergence  at  radius  r,  0^  Is  the  angu- 
lar divergence  without  thermal  effects,  and  (an/ST)  Is  the 
variation  of  the  Index  of  refraction  with  temperature. 

Letting  the  Index  n « 1,0  and  (1  - expC-az])  « az  (valid 
for  az  s 0.1),  the  equation  becomes 


, . e . _Pttz(an/aT)(i  - erf-arWl) 

r o 2Trkr 


(24) 


The  beam  Is  assumed  to  retain  a gausslan  distribution  as 
It  propagates.  In  order  to  use  a ray  transfer  matrix,  the 
beam  Is  now  assumed  to  be  an  Infinite  bundle  of  rays 
gausslanly  distributed  with  the  radius  r.  Each  ray  has  a 
given  divergence  dependent  on  r as  shown  by  equation  (24), 
with  the  beam  divergence  at  r = 0 and  r = «»  being  equal  to 
zero.  To  avoid  the  task  of  tracing  each  ray,  all  rays 
(except  at  r s 0)  are  now  assumed  to  have  some  average  diver- 


gence - 0Q.  Expanding  the  exponential  in  equation  (24) 
and  assuming  beam  truncation  by  the  aperture  at  the  exp(-2) 
point,  the  average  divergence  Is  calculated  to  be 


ft — Z-fT  - 0.0432Paz(an/ST) 

”r  ” ”o  ” ka 


(25) 
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for  no  wind. 

If  a wind  transverse  to  the  beam  path  Is  Introduced,  the 
steady-state  energy  equation  is  dlfflculy  to  solve  analytically. 
In  order  to  avoid  numerical  computer  rechnlques,  an  approximate 
solution  Is  used. 

As  the  beam  Is  assumed  to  have  a gaussian  Intensity  dis- 
tribution, It  Is  reasonable  to  assume  the  temperature  of  the 
beam  Is  also  gaussian  In  nature.  If  the  temperature  distri- 
bution with  X Is  assumed  to  be  of  the  form 


T = T exp{-2x^/a2)  = EXP 

o 


where  Is  the  temperature  at  x = 0,  then  the  first  and 
second  derivatives  are 


dx2  ar 


The  second  derivative  can  now  be  vrrltten  In  terras  of  the 
first  derivative  as 


^(1  4x. 

dx'x  - ^2^ 


Since  the  radius  of  Interest  Is  the  exp(-2)  point,  and  to 
simplify  matters  and  produce  an  equation  that  can  easily  be 
solved,  equation  (29)  Is  evaluated  at  the  point  x = a to  give 
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1(^ 

“ a'dx 


(30) 


i 

I 

! 

i 

i 


B^T 


dx 


2 - 


(Ref  7).  Then,  letting  (B^T/By^)  « (B^t/Bx'^), 
state  energy  equation  can  be  solved  to  give 

BT  ala 

dc^va  + 6k 
P 

Substituting  equation  (31)  Into  equatlon-.(22) , 
divergence  can  be  solved  as  before  to  be 


the  steady- 


(31) 


the  angular 


0. 1066  Paz(Bn/BT) 

a(dCpVa  + 6k) 


(32) 


for  the  angular  divergence  In  the  X direction. 

Any  wind  Is  assumed  to  reduce  blooming  only  In  the  vrlnd 

direction.  Neglecting  wind  steering  effects,  wind  reduces 

blooming  In  the  X direction  and  the  beam  becomes  elliptical 

In  shape  with  an  area  of  -rra  a , where  a and  a are  the  minor 

y * y 

and  major  radii  of  the  ellipse.  Total  beam  povrer  Is  calculated 
by  Integrating  over  the  beam  area.  To  simplify  the  Integration 
the  ellipse  Is  represented  by  a circle  of  the  same  area.  The 
equivalent  radius  of  the  circle  Is 

a = (aj^ay)^/^  ^ [(z"^)  (zS^)]^/^  (33) 

where  TZ  and  ” are  the  beam  divergences  In  the  X direction 

y 

and  the  Y direction,  respectively,  due  to  thermal  effects 

only.  Substituting  equations  (32)  and  (25)  for  ^ and  0” 

^ y 

gives 
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r ^ 


= 20 


/ 0.068  PazOn/aT)  \ 

' a(kdCpVa  + 6k^)^  / 


(34) 


where  the  absolute  value  of  (dn/3T)  Is  used. 
Thus  the  thermal  blooming  expression  is 


p _ 0.068  qzOn/ST) 

a(kdCpVa  + 


(35) 


and  the  transfer  matrix  for  calculating  beam  divergences 
only  is 


Complete  Transfer  Matrix 

Adding  matrices  (13)f  (20),  and  (36)  gives  -the  overall 
transfer  matrix  for  absorption,  turbulence,  and  thermal 
blooming.  The  final  form  is 


1 

r 

= 

N 

b( z)exp(-az) 
0 
0 


‘-I 


N-l 


N-1 


e 


N-l 


(37) 


and  is  the  form  that  should  be  used  for  all  calculations. 
If  a given  area  is  to  be  omitted,  for  example,  thermal 
blooming,  setting  6=0  would  give  the  transfer  matrix  for 
absorption  and  turbulence  only. 


Given  initial  values  of  P^,  x^,  and  0^,  equation  (37)  is 


applied  N times  (N  = L/z).  For  reasonable  results,  N should 
have  a value  greater  than  10,  and  is  restricted  by  the 


1 


PxT 

o ' 

o 

o 

J 

P.r  H 

1 

N 

N-l 

0 0 0 

^N-1 

(36) 

6 0 1 

_ N 

N-l 
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VALUES  FOR  d AND  (3n/3T) 


Altitude 

d X 10^ 

(an/ST)  X lO”"^ 

(km) 

(g/m^)* 

1 

o 

0.0 

1.22 

8.84 

0.5 

1.17 

8.52 

1.0 

1.11 

8.21 

1.5 

1.06 

7.90 

2.0 

1.01 

7.61 

2.5 

0.96 

7.32 

3.0 

0.91 

7.04 

4.0 

0.82 

6.50 

5.0 

0.74 

5.99 

8.0 

0.53 

4.63 

10.0 

0.4l 

3.85 

20.0 

0.09 

0.85 

*d  must  be  changed  to  units  of  g/crs?  when  used  In  equation 
(40) 


1 

) 

'i 


I 


earlier  assumption  that  az  s 0.1,  In  general,  az  £ 0,1  is 
valid  for  z ^ $00  m. 

After  applying  equation  (37)  N times,  the  final  beam 

power  P„,  final  radius  x„,  and  final  beam  divergence  e„  at 
N N N 

the  receiver  are  obtained.  To  get  a final  on-axls  intensity, 
it  is  recalled  that  the  total  power  in  a gaussian  beam  can 
be  calculated  by  integrating  over  the  volume  of  the  gaussian, 
or  from 


Zv  X 

Pjj  at  ^ ^ I^exp(-2r^/x^)  rdr  d0  (38) 

Where  is  the  on-axls  intensity  of  the  beam.  Equation 
(38)  can  be  solved  to  give 


I 


o 


N 


1.36  (x^) 


(39) 


which  is  used  to  determine  the  on-axls  intensity  at  the 
receiver. 

Concerning  the  final  transfer  matrix,  a few  comments  are 

necessary.  For  the  constants  in  the  matrix,  c^  and  k vary 

little  with  altitude,  and  the  values  c =1.0  J/g  °K  and  k = 

P 

0,024  J/m-sec®Kmay  be  used  at  all  times.  However,  the  values 
for  d and  (dn/3T)  decrease  with  altitude  and  should  be 
changed  for  every  kilometer  Increase  in  altitude.  Table  III 
is  a list  of  some  values  for  d and  (dn/3T).  The  values  for 
(3n/dT)  are  calculated  from 


hn  0.21d 
^ = T 


(40) 
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In  ®K"^  (Ref  81601)  where  d Is  In  g/cm^  and  T Is  in  °K  (Ref 

3«3-70). 
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Transmitter  rotation,  If  present,  Is  Included  In  the 
thermal  blooming  expression  by  solving  for  an  equivalent  wind 
velocity  of 

V ss  V •feus  (4l) 

o 

In  m/sec  where  v^  Is  the  transverse  wind  velocity,  Is  the 
transmitter  angular  velocity,  and  s Is  the  mean  distance 
from  the  transmitter  for  each  propagation  step. 


TABLE  IV 
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Comparison  of  Matrix  Code  and  COMBO 


Power 

Wind  (m/sec) 

(kw) 

Code 

1.0 

10.0 

25.0 

50.0 

Matrix 

I*  = 1. 14 

I a 6.07 

I a 12.17 

I a 20.89 

5 

X a 0,44 

X a 0.19 

X a 0.13 

X a 0.10 

COMBO 

I = 5.44 

X a 0,  12 

I a 14.01 

X a 0.12 

I a 15.12 

X a 0.12 

I = 15. 2R 
X a 0.  12 

Matrix 

I = 0.88 

I a 4.38 

I = 8.55 

I = 14.35 

10 

X a 0,70 

X a 0.32 

X a 0.22 

X a 0.17 

COMBO 

I a 6.29 

X a 0.26 

I a 24.09 

X a 0.13 

I = 28.73 

X a 0.12 

I a 30.24 
X a 0.  12 

Matrix 

I a 0.70 

I = 3.30 

I a 6.26 

I a 10.30 

20 

X a 1. 12 

X a 0.51 

X a 0.37 

X a 0.29 

COMBO 

I = 5.32 

X a 0.40 

I a 35.88 

X a 0.16 

r = 51.50 

X a 0.13 

I = 57.46 

X a 0.  12 

Matrix 

I a 0.54 

I a 2.38 

I a 4.40 

I a 7.06 

50 

X a 2.01 

X a 0.96 

X a 0.70 

X a 0.56 

j 

COMBO 

I = 3.03 

X a 0.85 

I a 54.42 

X a 0.20 

I = 89.71 

X a 0.  16 

I al20.47 
X a 0.13 

Matrix 

I a 0.45 

I = 1.93 

I = 3.50 

I = 5.53 

100 

X a 3.12 

X a 1.50 

X a 1.12 

X a 0.89 

COMBO 

I = 1.77 

X a 1.57 

I a 62.91 

X a 0,26 

I =122.86 

X a 0.19 

I =179.43 

X a 0.  l6 

*I  has  the  units  of  (10  ) w/m^ 
X has  the  units  of  m 


Ill,  Results  and  Discussion 


Comparison  with  COMBO 

Results  calculated  vising  the  matrix  code  were  originally 
intended  to  be  compared  with  real  data  collected  from 
controlled  experiments.  However,  such  data  Is  nearly  non- 
existent, especially  In  the  open  literature.  Therefore 
comparisons  were  made  with  the  computer  code  COMBO. 

A series  of  comparisons  were  made  using  various  Initial 
power  levels  and  wind  velocities.  The  results  of  the  trial 
cases  are  listed  In  Table  IV.  Except  for  the  cases  of  high 
powser  levels,  the  results  of  the  two  codes  are  consistent 
with  each  other. 

Comments 

In  Ollier  to  allow  a better  comparison  between  the  matrix 
code  and  COMBO,  a single  scenario  was  used  for  all  trial 
cases. 

The  problem  chosen  was  the  propagation  of  a laser  beam 
over  a path  5 km  long  at  a constant  altitude  of  6l  m.  The 
total  absorption  coefficient  was  the  same  for  both  codes. 

Since  COMBO  does  not  Include  scattering,  the  aerosol  scattering 
coefficient  b(z)  was  set  equal  to  1.0,  In  addition, 
turbulence  effects  were  not  included  In  the  trial  cases,  as 
the  increase  In  radius  due  to  turbulence  Is  small  (0.01  to 
2,0  cm),  and  the  two  codes  predicted  nearly  Identical  results 
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for  test  cases  Involving  turbulence  effects  only. 

The  difference  In  results  Is  thus  due  to  the  thermal 
blooming  models  of  the  two  codes.  For  lower  power  levels, 
results  predicted  by  the  two  codes  are  comparable.  For  high 
power  levels,  the  matrix  code  results  disagree  with  those  of 
COMBO.  The  assumption  that  mass  flow  In  the  X direction 
does  not  affect  blooming  In  the  Y direction  Is  no  longer 
valid. 

A more  exact  solution  of  the  energy  state  equation 
would  perhaps  correct  this  fault  In  the  matrix  model,  but 
the  obtalnment  of  that  solution  would,  most  likely.  Increase 
the  complexity  of  the  matrix  code  to  the  point  where  little 
would  be  gained  over  current  models  In  terms  of  computer 
time  and  storage. 

In  general,  the  major  area  of  disagreement  among  all 
models  Is  how  thermal  blooming  effects  are  calculated. 

Certainly  the  matrix  code,  if  It  does  not  predict  actual 
results,  at  least  calculates  the  worst  case. 

The  absorption  models  are  for  30°  N.  latitude.  To  adapt 
the  water  model  to  other  locations  requires  adjusting  the 
constants  A and  B In  equations  (9)  and  (10).  See  Appendix  A 
for  details.  If  absorption  coefficients  are  already 
available,  the  models  may  be  dropped  entirely. 

The  turbulence  model  predicted  the  same  results  as  that 
of  COMBO.  Two  cases,  one  at  61  m and  one  at  10,000  m,  were 
run  for  turbulence  only.  Both  the  matrix  code  and  COMBO 
calculated  a radius  Increase  of  1.8  cm  for  the  first  case, 
and  an  Increase  of  0,09  cm  for  the  second  case. 
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Although  all  the  examples  presented  assumed  a CO^  laser 
beam,  the  matrix  code  need  not  be  restricted  to  COg  lasers. 

By  properly  adjusting  (or  substituting)  the  absorption  coef- 
ficients, lasers  at  any  given  wavelength  may  be  traced  using 
the  matrix  formulation  developed  In  this  paper. 

Suggestions  and  Recommendations 

The  COg  model  Is  too  general  for  low  altitudes.  A more 
complete  model  based  on  local  conditions  Is  necessary  for 
best  results. 

The  thermal  blooming  model  appears  to  break  down  at  high 
power  levels.  A more  exact  solution  to  the  energy  state 
equation  can  be  Investigated  to  reduce  this  problem.  Another 
possibility,  one  less  likely  to  complicate  the  model,  is  the 
development  of  a correction  term  based  on  mass  flow  In  the  X 
direction  that  reduces  blooming  In  the  Y direction. 

The  model  Is  presently  restricted  to  a gausslan  beam. 

The  scope  of  the  model  would  be  greatly  Increased  If  a means 
of  Incorporating  other  beam  distributions  was  developed, 
possibly  by  representing  other  distributions  by  a sum  of 
gausslan  beams. 

There  Is  Interest  In  the  inverse  problem.  Given  a 
required  Intensity  at  the  receiver.  It  Is  desirable  to  know 
the  necessary  beam  power  at  the  transmitter.  If  the  matrix 
presented  here  can  be  solved  for  an  Inverse,  the  power 
needed  for  a specified  amount  of  energy  at  the  receiver 
could  be  calculated. 
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Appendix  A 
Water  Vapor  Model 

The  water  vapor  model  Is  derived  from  Ref  5.  The  model 
as  presented  Is  based  on  averages  and  Is  good  for  general 
results,  but  for  more  exact  results,  the  model  can  easily 
be  modified  to  reflect  local  conditions. 

The  authors  of  Ref  5 use  water  vapor  distributions  of 

WVP  = 11.0exp(-0.575H)  (42) 

for  January  and 

WVP  = 22.8exp(-0.505H)  (43) 

for  July  where  WVP  is  the  water  vapor  pressure  In  Torr. 

The  vertical  water  vapor  coefficients  B * 0.575  kni“^  and 
B = 0,505  km"^  are  from  plots  of  In(WVP)  versus  altitude  and 
are  the  slopes  of 
straight  lines  fitted 
to  data  (Ref  5 « 1475). 

See  Fig.  5.  If  local 
data  is  available,  a 
similar  plot  can  be 
mstde  and  new  values 
of  B calculated.  If 
not,  the  values  of  B 

given  earlier  may  be  Pig-  5.  Water  Vapor  Distribution 

Plot  (From  Ref  5 11475) 


^ 

used. 

The  constants  A t=  11.0  and  A = 22.8  Torr  are  based  on 
the  average  humidities  for  January  and  July.  However,  for 
any  given  temperature  and  relative  humidity,  A may  be  cal- 
culated directly  from 

A a 4. 58(RH)exp(5^23C“T5i^  - ^3)  (44) 

where  RH  Is  the  relative  hximldlty  (RH  s i.o)  and  T Is  the 
temperature  In  °K  (Ref  lil04). 


Appendix  B 

Polynomial  Curve  Pits  For  Fig.  2 

The  following  tables  contain  polynomials  that  are  the 

result  of  a curve  fitting  program  for  the  curves  In  Fig.  2. 

The  polynomials  are  used  to  calculate  the  CO^  absorption 

coefficient  a for  any  vertical  path  beginning  at  altitude 
c 

hj  and  ending  at  altitude  h^*  where  hj^  and  are  In  m,  by 
Integrating  the  polsrnoralal  (or  polynomials)  listed  for  the 
two  altitudes. 

The  polynomials  may  also  be  used  to  calculate  a for 

c 

any  given  altitude  by  substituting  the  altitude  h In  m 
directly  Into  the  correct  polynomial. 
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TABLE  VI 

July  Polynomial  Curve  Fit  For  Fig.  2 


Altitude 

Polynomial 

(km) 

0-15  8.535241  X 10“^  - 1.439371  r,  10“^h  + 2.2556OI 

X 10” V - 3.032077  X 10“^^^  + 2.079833 
X 10“ - 5.087308  X i0‘^^h5 

15.47  4,409444  X 10-2  „ 8.501794  :c  10“^h  + 6.403721 

X 10"^ V - 2.  137123  X 10" + 3.280703 
X 10” - 1.912325  X 10"2^V.5 

47-51  -5.381483  X 10"^  + 3.443774  X 10"^h  - 7.221806 

X 10"^®h^  + 4.994872  x 10”^^h^ 

51-65  8.595239  X 10"^  - 5,584211  X 10"^h  + 1.371591 

X 10”^h^  - 1.5063  X 10" + 6.229936  X 10"^°h^  . 


» 


TABLE  V 

January  Poljmomlal  Curve  Fit  For  Fig.  2 


Altitude 

(km) 

Polynomial 

0—2 

6.177301  x 10“^  - 4.243258  X 10“^h  + 3.331791 

^ X 10“^\i^ 

2—12 

7.1772259  x 10“^  - 1. 018966  X 10“^h  + 5.438857 

X 10“^^h^  - 1.440797 X 10“ + 2. 111184 x 10" 

12—17 

2.168201x  10“^  - 1,708914  X 10“^h  + 4.794412 

X 10“^^h^  - 4.821832 X 10"^^h^ 

17—4? 

2.749'+15x  10“^  - 6.043183  X 10“^h  + 4.952857 

X 10“^%^  - 1.73234  X 10“ + 2.748747  X 10“^ 

- 1.644807  X 10“^'’'h^ 

^7—51 

-5.0246274x  10“^  + 3.209868  X i0“^h  - 6.729368 

X 10“^°h^  + 4.656468 X 10“^^h^ 

51—65 

1.2415226x  10“^  - 5.636995  X 10“^h  + 8.609204 

X 10“^^h^  - 4.4187473 X 10“^^h^ 

h is  the  altitude  in  meters  (From  Ref  9 >1552) 
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Sample  Calculator  Programs 


The  following  pages  contain  two  sample  calculator  programs 
for  a single  propagation  scenario  as  an  example  of  how  the 
matrix  formulation  may  be  used  on  a hand  calculator.  The 
first  program  Is  designed  for  a Texas  Instruments  SR  ^6,. 

The  second  program  Is  designed  for  a Hewlett  Packard  HP  25. 

The  propagation  scenario  that  the  programs  solve  Is  the 
coaltitude  transmission  of  a laser  beam  from  a transmitter 
to  a receiver  over  a distance  L In  propagation  steps  of 
distance  z.  The  beam  has  some  Initial  power  P and  radius  x 
and  Is  assimed  t'o  be  focused  at  a point  distance  L away. 

Because  the  effect  of  turbulence  Is  small  compared  to 
thermal  bDoomlngt  It  Is  not  Included  In  these  programs.  If 
desired,  however,  It  may  be  added. 
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SR  56  Program 

The  memory  locations  have  the  following  values  1 

(0)  N — number  of  propagation  steps 

(1)  P*””  (0. 865) (Initial  beam  power)  due  to  aperture 

truncation  at  the  exp(-2)  point,  In  watts 

(2)  X — Initial  beam  radius  In  meters 

o 

(3)  Initial  beam  divergence  In  radians 

(4)  -a~  negative  value  of  total  absorption  coefficient,  km' 

(5)  3C  --  Initial  beam  radius  In  meters 

o 

(6)  Cj^—  constant  = 0.  OOO68  z(dn/aT) 

(7)  Cg--  constant  = 0,024vd 

(8)  C^—  constant  = 6k^  = O.OO3456 

(9)  z --  propagation  step  distance  In  meters 


The  program  keystrokes  arei 
LRN 

RCL  3 X RCL  9 + RCL  2 = STO  2 

RCL  3 - RCL  1 X RCL  4 x RCL  6 
4-((RCL  7 X RCL  5 + RCL  8)v5r 
X RCL  5)  = STO  3 

RCL  2 
STO  5 

RCL  1 X RCL  4 e*  a STO  1 
2nd  dsz  00 

RCL  1 (1.36  X RCL  2 y?)  a 

R/S 

RST 

LRN 


Calculates  new 
radius 

Calculates  new 
beam  divergence 


Replaces  old  value 
of  radius  with  new 

Calculates  new 
power 

Do  loop  Instruction 

Calculates  final 
on>axls  Intensity 


Upon  pressing  the  R/S  key  the  program  will  run  to  com- 
pletion, displaying  the  on-axls  Intensity  at  the  receiver. 
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The  final  beam  power  Is  stored  In  memory  location  (1)  and 
the  final  beam  radius  is  stored  In  memory  location  (2). 


Example  I For  N a 50,  L = 5000  m,  z = 100  m,  P = $000  w, 

h a 61  m,  a 0.2  m,  a a 0.07562  km“^,  (dn/dT)  a S.S^xlO""^, 

V a 1.0  m/sec,  d a 1.22 x 10^  g/m^,  the  final  values  are 

Pjj  a 2960  w,  Xjj  a 0.47  m,  and  I a ll400  w/m^. 

The  initial  values  are  P^  a 4325  w,  x^  a 0.2  m,  and 

6^  a >4. Ox  10“^  radians, 
o 
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HP  25  Program 

The  memory  locations  have  the  following  values t 

(0)  az—  total  absorption  coefficient  In  lan“^xz  In  km 

(1)  (0« 865) (Initial  beam  power)  due  to  aperture 
® truncation  at  the  ezp(-2)  point 

(2)  z — Initial  beam  radius  In  meters 

o 

(3)  0 — Initial  beam  divergence  In  radians 

o 

(4)  0 — Initial  value  for  do  loop 

(5)  Cj^—  constant  = 0. OOO68  z(an/aT) 

(6)  Cg—  constant  Cg  = 0.024vd 

(7)  constant  =*  6k^  s 0.003^56 


The  program  keystrokes  arei 

RCL  1 ( 

RCL  0 i 

z 

RCL  5 
z 

RCL  6 
RCL  2 
z 

RCL  7 
+ 

nT* 

RCL  2 

z 


Calculates  new  divergence 
angle 


RCt  3 


RCL  3 

STO  3 

CLR  X 
1 
0 
0 
z 

RCL  2 
+ 

STO  2 


Stores  new  divergence  angle 


Calculates  new  radius  1 In 
place  of  100,  enter  step 
distance  z 
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Calculates  new  power 


RCL  0 
CHS 
e* 

RCL  1 

X 

STO  1 

1 Do  loop  Instructions  I In  place 

STO  +4  of  50t  enter  the  number  of 

RCL  4 propagation  steps 

5 

0 

X Y 
GO  TO  01 

GO  TO  00  Ends  program 

Upon  pressing  the  R/S  key  the  program  will  run  to  com- 

pletion, storing  the  final  power  In  memory  (1)  and  the 
final  radius  in  memory  (2),  To  calculate  the  final  on-axls 
Intensity,  the  keystrokes  are 

RCL  1 
RCL  2 

t 

1.36 


The  display  will  show  the  final  Intensity. 
Example  I See  example  for  SR  56. 
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